0.1 INTRODUCTION

This analysis assesses the slopes for injury scores over time between placebo and felzartamab arms. We used a linear mixed-model to account for repeated measurements over time.

0.2 DATA WRANGLING


First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# CRAN packages
library(tidyverse)
library(lme4)
library(lmerTest)
library(performance)
library(ggeffects)
library(flextable)
library(ggpubr) 
library(patchwork) 
library(ggprism) 
# Bioconductor libraries
library(Biobase) # BiocManager::install("Biobase")
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
source("C:/R/CD38-effect-of-treatment/natmed/code/functions/get_slope_function.r")
# load reference set
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
# load refression plots
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_regression.RData")
# load violin plots
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")

# Seed for reproducibility
set.seed(42)


# DEFINE THE DATA ####
data <- set %>%
    pData() %>%
    dplyr::select(Center, Patient, Felzartamab, Group, IRRAT30, IRITD3, IRITD5) %>%
    dplyr::filter(Patient %nin% c(15, 18)) %>%
    mutate(
        time = case_when(
            Group == "Index" ~ 0,
            Group == "FU1" ~ 0.46,
            Group == "FU2" ~ 0.99
        ),
        linetype = "solid" %>% factor()
    )


Now we fit the linear mixed-models.

# FIT MODELS ####
fit_IRRAT30 <- lmer(IRRAT30 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD3 <- lmer(IRITD3 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD5 <- lmer(IRITD5 ~ Felzartamab * time + (time | Patient), data)


We check the validatity of the model assumptions for IRRAT30.

# test model assumptions
performance::check_model(fit_IRRAT30)

performance::check_normality(fit_IRRAT30)
## OK: residuals appear as normally distributed (p = 0.605).


We check the validatity of the model assumptions for IRITD3.

# test model assumptions
performance::check_model(fit_IRITD3)

performance::check_normality(fit_IRITD3)
## OK: residuals appear as normally distributed (p = 0.070).


We check the validatity of the model assumptions for IRITD5.

# test model assumptions
performance::check_model(fit_IRITD5)

performance::check_normality(fit_IRITD5)
## OK: residuals appear as normally distributed (p = 0.090).

0.3 ORGANIZE THE RESULTS


The next step is to extract and summarize the model results.

# EXTRACT THE SLOPES FROM THE MODEL FITS ####
# IRRAT30
# Create list object into which to place model results
slopes_IRRAT30 <- list()
# Acute slopes
slopes_IRRAT30[[1]] <- get_slope(
    .model_obj = fit_IRRAT30, .name = "Placebo",
    .contrasts = c(0, 0, 1, 0) #
)
slopes_IRRAT30[[2]] <- get_slope(
    .model_obj = fit_IRRAT30, .name = "Felzartamab",
    .contrasts = c(0, 0, 1, 1)
)
slopes_IRRAT30[[3]] <- get_slope(
    .model_obj = fit_IRRAT30, .name = "Intergroup Difference",
    .contrasts = c(0, 0, 0, 1)
)
model_IRRAT30 <- slopes_IRRAT30 %>%
    bind_rows() %>%
    mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
    dplyr::select(name, result, p_value, p_adjusted)

# IRITD3
# Create list object into which to place model results
slopes_IRITD3 <- list()
# Acute slopes
slopes_IRITD3[[1]] <- get_slope(
    .model_obj = fit_IRITD3, .name = "Placebo",
    .contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD3[[2]] <- get_slope(
    .model_obj = fit_IRITD3, .name = "Felzartamab",
    .contrasts = c(0, 0, 1, 1)
)
slopes_IRITD3[[3]] <- get_slope(
    .model_obj = fit_IRITD3, .name = "Intergroup Difference",
    .contrasts = c(0, 0, 0, 1)
)
model_IRITD3 <- slopes_IRITD3 %>%
    bind_rows() %>%
    mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
    dplyr::select(name, result, p_value, p_adjusted)

# IRITD5
# Create list object into which to place model results
slopes_IRITD5 <- list()
# Acute slopes
slopes_IRITD5[[1]] <- get_slope(
    .model_obj = fit_IRITD5, .name = "Placebo",
    .contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD5[[2]] <- get_slope(
    .model_obj = fit_IRITD5, .name = "Felzartamab",
    .contrasts = c(0, 0, 1, 1)
)
slopes_IRITD5[[3]] <- get_slope(
    .model_obj = fit_IRITD5, .name = "Intergroup Difference",
    .contrasts = c(0, 0, 0, 1)
)
model_IRITD5 <- slopes_IRITD5 %>%
    bind_rows() %>%
    mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
    dplyr::select(name, result, p_value, p_adjusted)

0.4 SUMMARIZE RESULTS

1 CREATE TABLE SUMMARIES OF SLOPES


IRRAT30.

size_font <- 10
model_IRRAT30 %>%
    dplyr::select(-p_value) %>%
    flextable::qflextable() %>%
    flextable::set_table_properties(layout = "autofit") %>%
    flextable::set_header_labels(
        name = "Group", result = "IRRAT30 Slope (95%CI)",
        "p_adjusted" = "FDR"
    ) %>%
    flextable::bold(i = NULL, part = "header") %>%
    flextable::fontsize(size = size_font, part = "all")

Group

IRRAT30 Slope (95%CI)

FDR

Placebo

0.32 (-0.04 to 0.67)

0.1133

Felzartamab

-0.29 (-0.64 to 0.07)

0.1133

Intergroup Difference

-0.60 (-1.10 to -0.10)

0.0546


IRITD3.

model_IRITD3 %>%
    dplyr::select(-p_value) %>%
    flextable::qflextable() %>%
    flextable::set_table_properties(layout = "autofit") %>%
    flextable::set_header_labels(
        name = "Group", result = "IRITD3 Slope (95%CI)",
        "p_adjusted" = "FDR"
    ) %>%
    flextable::bold(i = NULL, part = "header") %>%
    flextable::fontsize(size = size_font, part = "all")

Group

IRITD3 Slope (95%CI)

FDR

Placebo

0.06 (-0.06 to 0.18)

0.3453

Felzartamab

-0.14 (-0.26 to -0.02)

0.0408

Intergroup Difference

-0.20 (-0.37 to -0.02)

0.0408


IRITD5.

model_IRITD5 %>%
    dplyr::select(-p_value) %>%
    flextable::qflextable() %>%
    flextable::set_table_properties(layout = "autofit") %>%
    flextable::set_header_labels(
        name = "Group", result = "IRITD5 Slope (95%CI)",
        "p_adjusted" = "FDR"
    ) %>%
    flextable::bold(i = NULL, part = "header") %>%
    flextable::fontsize(size = size_font, part = "all")

Group

IRITD5 Slope (95%CI)

FDR

Placebo

0.09 (-0.04 to 0.22)

0.15720

Felzartamab

-0.13 (-0.26 to -0.01)

0.05895

Intergroup Difference

-0.23 (-0.41 to -0.05)

0.04200

1.1 PLOT THE RESULTS

# DEFINE LEGEND FOR VIOLIN PLOTS ####
panel_legend_violin <- plots_artANOVA %>%
    dplyr::filter(variable == "AMAT1") %>%
    pull(plot_violin) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot() +
    theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))


# DEFINE LEGEND FOR REGRESSION PLOTS ####
panel_legend_regression <- plots_regression %>%
    dplyr::filter(variable == "IRRAT30") %>%
    pull(plot) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot()


# DEFIN JOINT LEGEND FOR REGRESSION AND VIOLIN PLOTS ####
panel_legend_all_injury <- panel_legend_violin + panel_legend_regression


# DEFINE INJURY VIOLIN PLOTS ####
plots_violin <- plots_artANOVA %>%
    dplyr::filter(category %in% c("injury")) %>%
    pull(plot_violin)


# MAKE PANEL OF SLOPE AND VIOLIN PLOTS ####
plot_panels_all_injury <- patchwork::wrap_plots(
    plots_violin[[1]], plots_violin[[2]], plots_violin[[3]],
    plots_regression$plot[[1]], plots_regression$plot[[2]], plots_regression$plot[[3]],
    plots_regression$grob_table[[1]], plots_regression$grob_table[[2]], plots_regression$grob_table[[3]],
    nrow = 3,
    ncol = 3,
    guides = "collect",
    heights = c(1, 1, 0.5)
) +
    patchwork::plot_annotation(tag_levels = list(c(LETTERS[1:6], rep("", 3)))) &
    theme(
        legend.position = "none",
        # plot.title = element_blank(),
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        plot.tag = element_text(size = 20, face = "bold", vjust = 1)
    )


ggarrange(
    panel_legend_all_injury,
    plot_panels_all_injury,
    nrow = 2,
    heights = c(0.15, 1.5)
)